#include "mrate.h"
#include <cmath>
#include "ext_rational.h"

const double f1_gamma(const double x, const double rlen){
  return x + 1.0 / pow(x, rlen-1) / (x-1);
}

const double f2_gamma(const double x, const double rlen){
  return x + 1.0 / pow(x, 2*rlen-1) / (x*x-1);
}

const double calc_mrate(const input_data & INPUT, const grate_config & CONFIG, const double grate){
  double res = (CONFIG.PRECISE_LOWER) ? search_mrate(INPUT, grate, &f2_gamma) : search_mrate(INPUT, grate, f1_gamma);
  return (res == 0) ? -2 : res;
}

const double search_mrate(const input_data & INPUT, const double grate, FGAMMA f){
  double rlen = INPUT.ROOT_LENGTH;
  // ternary search to calculate f(gamma)
  double l = 1, r = grate, p, q;
  while (r - l > INPUT.delta){
    p = (2*l+r)/3.0;
    q = (l+2*r)/3.0;
    if (f(p, rlen) < f(q, rlen)) r = q;
    else l = p;
  }
  l = (l+r)/2.0;
  r = grate;
  // Now l is the minimum of f(gamma)
  // binary search to find the place where f(gamma) = grate
  if (f(l, rlen) > grate){
    return 0;
  }
  while (r - l > INPUT.delta){
    q = (l+r)/2.0;
    if (f(q, rlen) < grate) l = q;
    else r = q;
  }
  return (l+r)/2.0;
}
